-9- 



REMARKS 

The present Amendment makes editorial changes in the title, 
specification and claims, and adds an Abstract to conform the present PCT 
application to the requirements of United States patent practice. The 
5 cancellation of claims 1-13 in favor of the claims presented herein has been 
done solely because the amount of bracketing and underlining that would 
have been necessary to conform those claims to the requirements of United 
States patent practice would have been unduly burdensome and confusing. 
No change in the claim language has been made for the purpose of 
10 distinguishing any of the claims over the teachings of the prior art of record. 
Accordingly, no change in the claim language is considered by the Applicants 
as a surrender of any of the subject matter encompassed within the scope of 
the original claims 1-4. 

Early consideration on the merits is respectfully requested. 
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APPARATUS FOR ANALYS I NG CARDIAC EVENTS 
T e chn i cal fi e ld 

SPECIFICATION 
TITLE 

5 "APPARATUS FOR ANALYZING CARDIAC EVENTS" 

BACKGROUND OF THE INVENTION 

Field of the Invention 

The present invention relates to an apparatus for analys i ng cardiac 
events detected in electrograms, EGMs, and to a heart stimulator provided 

10 with such an apparatus. 

In the following the expression card i ac evont "cardiac event" denotes 
the depolarization phase in the cardiac cycle^ which for atrial signals is 
commonly known as P wave and for ventricular signals as R wave or QRS 
complex. 

15 Background 

Description of the Prior Art 

In the field of devices for cardiac rhythm management (CRM), accurate 
rhythm classification is an increasingly important aspect. Pacemakers are 
primarily used to assist in bradycardia or when the electrical propagation path 

20 is blocked, whereas the primary use of implantable cardioverter defibrillators 
(ICD) is to terminate ventricular arrhythmia, a life-threatening condition if not 
acute l y immediatelv treated. In both types of devices, accurate event 
classification of the electrogram signal is needed for identifying, e.g., atrial 
and ventricular fibrillation in order to give appropriate therapy for the detected 

25 arrhythmia. For pacemakers, this may imply necessitate changing the pacing 
mode in order to stabilize the ventricular rhythm during an episode of atrial 
fibrillation. An ICD responds to ventricular fibrillation by giving a dofibri ll atory 
defibrillatinq shock intended to terminate the fibrillation. 

Ever increasing demands are put on both kinds of devices to better 

30 handle their primary task as well as to manage other tasks than those 
originally intended for. One such task may be, for an implantable medical 
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device, be to identify atrial flutter in order to terminate it by atrial pacing or to 
defibrillate atrial fibrillation. Although it is not a life-threatening arrhythmia, 
atrial fibrillation is an inconvenience to the patient and Increases the risk for 
other diseases such as stroke. Atrial pacing may also be one way of 
5 terminating supraventricular tachycardias. An ICD specific task is to identify 
atrial fibrillation in order to not mistake it for ventricular fibrillation and the risk 
of giving an unnecessary, and possibly harmful, defibrillation shock. Another, 
more general, utilization is to efficiently store rhythm data for later analysis 
and evaluation, already done in modern ICD's. By collecting data, better 

10 knowledge of the evolution of cardiac diseases and the functionality of the 
device can be obtained. 

Clustering represents an important task within the classification 
problem where each individual event is assigned to a cluster of events with 
similar features. Labelling of the clusters, i.e., associating the cluster with a 

15 specific cardiac rhythm, completes the classification such that the device can 
provide proper therapy when needed. However, certain constraints 
distinguish clustering in CRM devices from clustering in general. In order to 
give immediate therapy, it requires clustering to be done in real-time, thus 
excluding many iterative clustering algorithms such as k-means clustering and 

20 competitive learning. Various methods have recently been presented 
concerning clustering of signals from the surface electrocardiogram (ECG), 
based on, e.g., self-organizing maps or fuzzy hybrid neural networks, see M. 
Lagerholm, C. Peterson, G. Braccini, L. Edenbrabdt, and L. Sornmo, 
Clustering ECG complexes using Hermite functions and self-organizing 

25 maps", IEEE Trans. Biomed. Eng., vol. 47, pp. 838-848, July 2000, and S. 
Osowski and T. Linh, "ECG beat recognition using fuzzy hybrid neural 
network", IEEE Trans. Biomed. Eng., vol. 48, pp. 1265 -1271, November 
2001. However, most clustering algorithms used for ECG analysis are 
computationally rather complex and therefore unsuitable for implantable CRM 

30 devices. Furthermore, not much a priori morphologic information is 
associated with the various rhythms in the electrogram (EGM); this is in 
contrast to the more well-defined ECG. 



■ t 

-3- 

MARKED-UP COPY 

Pr e v i ous l y pr e s e nt e d Previous work In the area of intracardiac event 
classification mainly focus focused on discrimination of a specific condition in 
order to discern, e.g., atrial fibrillation from other atrial tachyarrythmias, see 
A. Schoenwald, A. Sahakian, and S. Swiryn, "Discrimination of atrial fibrillation 
5 from regular atrial rhythms by spatial precision of local activation direction", 
IEEE Trans, Biomed. Eng., vol. 44, pp. 958-963. October 1997. Other 
applications involve discrimination of ventricular from supraventricular 
tachycardia, see L. Koyrakh, J. Gillberg, and N. Wood, "Wavelet based 
algorithms for EGM morphology discrimination for implantable ICDs", in Proc. 

10 Of Comp. In Card. (Piscataway, NJ, USA), pp. 343 - 346, IEEE, IEEE Press. 
1999, and G. Gronefeld, B. Schulte. S. Hohnloser, H. - J. Trappe, T. Korte, 
C. Stellbrink, W. Jung, M. Meesmann, D. Bocker, D. Grosse - Meininghaus, 
J. Vogt, and J, Neuzner, "Morphology discrimination: A beat-to-beat algorithm 
for the discrimination of ventricular from supraventricular tachycardia by 

15 implantable cardioverter defibrillators", J. Pacing Clin. Electrophysiol., vol. 24, 
pp. 1519 - 1524, October 2001 . More general classification algorithms, which 
in turn involve training on individual patients, have been based on ana l ogu e 
analog neural networks or wavelet analysis for morphologic discrimination of 
arrhythmias. 

20 PCT Application WO 97 39 681 describes a defibrillator control system 

comprising a pattern recognition system. The intracardiac electrogram signal 
is digitised and delivered for feature selection into a selector. The feature 
selector outputs selected features to a trained classifier to provide information 
as to what group the produced signal should be clustered, e.g. ventricular 

25 tachycardia. The classifier outputs the classified information for use for a 
therapeutic decision. 

US 5 271 ^11 United States Patent No. 5,271,411 discloses an ECG 
signal analysis and cardiac arrhythmia detection by extraction of features from 
a scalar signal. A QRS pattern vector is then transformed into features 

30 describing the QRS morphology, viz. a QRS feature vector. A normal QRS 
complex is identified based on the population of QRS complexes located 
within clusters of QRS features within a feature space having a number of 
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dimensions equal to the number of extracted features. The extracted 
morphology information is then used forjudging whether a heartbeat is normal 
or abnormal, 

US 5 638 823 United States Patent No. 5,638,823 describes non- 
5 invasively detecting of coronary artery diseases. A wavelet transform is 
performed on an acoustic signal representing one or more sound event 
caused by turbulence of blood flowing in an artery to provide parameters for a 
feature vector. This feature vector is used as one input to neural networks, 
the outputs of which represent a diagnosis of coronary stenosis in a patient. 

10 In Michael A. Unser et al, "Wavelet Applications in Signal and Image 

Processing IV", Proceedings SPIE - The International Society for Optical 
Engineering, 6-9 August 1996. vol. 2825, part two of two parts, pp. 812 - 821, 
a wavelet packet based compression scheme for single lead EGGs is 
disclosed, including QRS clustering and grouping of heart beats of similar 

15 structures. For each heart beat detected, its QRS complex is compared to 
templates of previously established groups. Point-by-point differences are 
used as similarity measures. The current beat is assigned to the group whose 
template is most similar, provided predetermined conditions are satisfied. 
OthenA^ise a new group is created with the current QRS complex used as the 

20 initial group template. 

SUMMARY OF THE INVENTION 

Futur e pac e mak e rs and ICD's wil l includo more advanced moons for 
arrhythmia d e t e ct i on and therapy, and th e purpos e An obiect of the present 
invention is to propos e provide a technique for separation of cardiac rhythms 
25 in a reliable way on the basis of electrogram event clustering 
D i sc l osur e of th e i nv e nt i on 

The above purpose i s obtain e d by an apparatus according to c l aim 1 . 
The above obiect is achieved in accordance with the invention by an 
apparatus for analyzing cardiac events detected in electrogram (EGMs) 
30 having a feature extraction unit that derives features of the cardiac events for 
discriminating among different tvpes of detected cardiac events, and a 
clustering unit that groups cardiac events with similar features into a cluster. 
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defined by predetermined cluster features, wherein the feature extraction unit 
determines a feature vector that describes waveform characteristics of the 
cardiac event EGM signals bv means of a wavelet transform, and wherein the 
clustering unit determines a distance or distances between the feature vector 
5 and corresponding cluster feature vectors in order to assign the cardiac event 
in guestion to the cluster that results in a minimum distance, as long as the 
minimum distance is less than a predetermined threshold value. 

Certain arrhythmias are diagnosed immediately to give proper therapy, 
while, for others, it may be sufficient to record the rhythm for data collection 

10 purposes. Thus, the classification problem, viz. to label the rhythms based on 
clusters using clinical terms, may not always be necessary to implement. 

According to advantag e ous e mbodiments In an embodiment of the 
apparatus according to the invention both morphologic and temporal data are 
considered for clustering. Morphologic features are efficiently extracted by 

15 use of the dyadic wavelet transform after which the events are grouped by a 
leader-follower clustering embodiment. The event detection problem, based 
on the same transform, is previously treated in Swedish patent application no. 
0103562-5. 

According to another advantageous embodiment of the apparatus 
20 according to the invention an i nt e grat i ng m e ans integrator is provided to 
integrate said distance over a predetermined period of time. By integrating 
the distance over a period of time it is possible to distinguish irregular 
rhythms, like atria! fibrillation, from regular rhythms. The integral total 
distance in the case of atrial fibrillation will be high whereas regular rhythms 
25 will result in a lower total distance. 

The invention also relates to a heart stimulator provided with the 
above-mentioned apparatus for controlling the therapeutic stimulation 
depending on arrhythmia detection. 
Brief description oftho drawings 
30 To e xpla i n th e inv e nt i on in gr e at e r d e tai l e mbod i m e nts of th e inv e nt i on, 

chos e n as e xampl e s, w i l l now b e described w i th roforonco to the onc l os e d 
drawings, on which figure 1 shows impu l se r e spons e s of a f il t e r bank us e d for 
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cardiac ovont dotoction, figure 2 i s a flow chart of clust e ring a l gorithm p e rform 
with tho apparatus accord i ng to tho inv e ntion, figure 3 ill ustrates th e 
computational complex i ty for d i ffer e nt clust e ring a l gor i thms, figure A pr e sents 
in (a) clust e r i ng p e rformanc e in torms of probability of a corr e ct clust e r i ng of 
5 an e v e nt, Pcg and probabi li ty of a dominant e v e nt i n a clust e r Pos and i n (b) 
th e spr e ad in th e numb e r of clusters, figur e 5 i l l ustrat e s i n (a) clust e ring 
p e rformanc e in terms of Pqe and Pgg as a function of a distanc e threshold q 
and in (b) th e numb e rs of clust e rs as a funct i on of n, figur e 6 e x e mplifi e s th e 
c l ust e ring p e rformance for one set of param e t e rs, figur e 7 e x e mp li f ie s 

10 clust e ring r e sult on a concat e nat e d EGM for throe cases, figure 8 is a flow 
chart i l lustrating th e funct i on in broad out l ine of an ombod i mont of th e 
apparatus accord i ng to th e inv e ntion, and figur e 9 is a b l ock diagram of an 
exemplifying ombod i mont of a h e art stimu l ator provided w i th an apparatus 
according to th e i nv e ntion. 

15 Description of preferred e mbodim e nts 

DESCRIPTION OF THE DRAWINGS 
Figures la and lb illustrate impulse responses of a filter bank used for 

cardiac event detection in accordance with the present invention. 

Figure 2 is a flow chart of a clustering algorithm performed in the 

20 apparatus according to the present invention. 

Figures 3a and 3b illustrate the computational complexity for different 

clustering algorithms. 

Figure 4a illustrates clustering performance in terms of probability for 

correct clustering of an event, as well as the probability of a dominant event in 

25 the cluster and Figure 4b illustrates a histogram spread of the number of 

clusters. 
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Figure 5a illustrates clustering performance as a function of a distance 
threshold, and Figure 5b illustrates the number of clusters as a function of the 
distance threshold. 

Figure 6 illustrates the clustering performance for one set of 
5 parameters- 
Figures 7a and 7b are examples of the clustering result for a connected 
EGM for three cases. 

Figure 8 is a flowchart describing the general functioning of an 
embodiment of an apparatus according to the invention, 
10 Figure 9 is a block diagram of a heart stimulator providing with a 

cardiac event detecting apparatus in accordance with the present invention. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 

P wave detection and feature extraction 

A signal model assuming that the event waveform is composed of a 
linear combination of representative signals is considered. The feature 
extraction problem is then to estimate the individual components of the 
representative signals since each morphology will have its own linear 
combination. By using the dyadic wavelet transform, different widths of the 
two fundamental monophasic and biphasic waveforms are included in the 
model at a low cost. 
Feature extraction 

It is assumed that the QRS waveform is composed of a linear 
combination of representative signals, 

H = [hi hp] (2) 

where each function, hj, j= 1. P, is of size M x 1. The only restriction on 

H is that it must have full rank. Different morphologies, s(n), are modelled by 
the P X 1 coefficient vector 9(n), with the linear model 

s(n) =He(n) (3) 



15 



20 



25 
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where n is a temporal variable describing when the event occurs. The 
observed signal. x(n), is assumed to be modelled by, 

x(n) =He(n) + w(n) (4) 
where the M x 1 noise vector w(n) is assumed to be zero-mean, white, and 

2 

Gaussian with variance a — . Consequently, x(n) is defined as 

CO 



x(n) = 



x(n) 



(5) 



x(n + M-l) 

implying an event duration of M samples, beginning at n. The probability 
density function x(n) for a specific realization of 0(n), p(x(n); 0(n)), is thus 
given by 



p(x(n)0(n))= 



1 



2Tra 



2 IM 



exp 



1 



-(x(n)-He(n)r(x(n)-He(n)) 



2a 



CO 



(6) 



In this model, a complete description of a QRS complex is provided by 
the deterministic unknown parameter vector 0(n). The absence of a QRS 
complex corresponds to the case when 0(n) is equal to 0, where 0 is the P x 1 
zero vector. In general, no a priori knowledge is available on 6(n), and 
therefore an estimate is required before detection can take place. 
Furthermore, only one event is assumed to take place within the observation 
interval 0 < n < N. 
Fiiterbank representation 

The descriptive functions in H have been selected such that the 
following three aspects have been taken into particular account: 

1. the main morphologies of the QRS complex are mono- and/or 
biphasic 

2. the broad range of QRS complex durations, and 

3. a low complexity implementation. 

The wavelet transform is particularly suitable since it is a local 
transform, i.e., it provides information about the local behaviour of a signal. 
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10 



One wavelet decomposition method which may be efficiently implemented is 
the dyadic wavelet transform. By careful selection of the filters, a suitable 
filter bank including mono-and biphasic Impulse responses can be obtained. 
A symmetric lowpass filter, f(n). is used repeatedly in order to achieve proper 
frequency bands. This filter is combined with one of two filters, gb(n) or gm(n), 
which together define the waveform morphology (the subindices b and m 
denote bi- and monophasic, respectively). For the biphasic case, the 
recursion is expressed as, 

hi.b(n)= gb(n) 

h2.b(n)= f(n) * gb(2n) 

h3.b(n)= f(n)*f(2n)*gb(4n) 



= f(n)* ... * 1(2^"- ^n)* g,(2^-^n ) 



(7) 



in which the subindex qmax represents the maximum (coarsest) scale. 
15 It is now possible to present an expression for-H which is composed of 

one biphasic and one monophasic part. 



H=Lh. H.J 



(8) 



20 



where the biphasic Hb is defined as 



^(M-l)... ^(M-1) 

min'" max'" 



(9) 



where the subindex qmin represents the minimum scale and qmin < qmax. The 
25 monophasic matrix Hm is computed in a corresponding way. The reversed 
order of the columns in H , denoted with H in (8), is introduced in order to be 
consistent with the model assumed in (4). 

In order to mimic the desired mono- and biphasic waveforms with a low 
complexity filter bank structure, short filters with small integer coefficients 
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were used. In (7), the impulse response f(n) was chosen as a third order 
binomial function, 

F(z) = (l + 1 + 3z-' + 3z-' + Z-' (1 0) 

where F(z) is the Z-transform of f(n). For the biphasic filter bank, the 
5 filter gb(n) was selected as the first order difference, 

gb(") = [-1 1] 01) 
The filter gm(n) was chosen such that a compromise between the 
requirement of a DC gain equal to zero and an approximately monophasic 
impulse response was achieved. A reduction of complexity results when gb 
10 (n) is reused, 

g^(n) = g,(n)*g,(n)=[l-2l] (12) 

For this particular choice of gm(n) and by using Mallafs algorithm, it is 
possible to calculate both the biphasic and the monophasic filter output from 
each scale by using f(n) once and gb(n) twice. The filter bank impulse 

15 responses are shown in Fig. 1 Figs, la and 1b . The filter bank includes two 
orthogonal signal sets where the width of the signal varies within each set. In 

(a) Fig, la . hj,b (n) is shown for j = 2, ,4 from top to bottom, and in (b) Fig, 

lb, the corresponding hj.m(n) are shown. 
ML parameter estimation 

20 The unknown coefficient vector e(n) can be estimated by using the 

maximum likelihood criterion according to, 

0(n) = arg_p(x(n)e(n)) (13) 
e(n) 

However, e(n) is only of interest for those n for which the probability of 
an event, or, equivalently, for which the likelihood ratio test function T(x(n)) is 
25 maximized, 

n = arg,3,T(x(n)) (14) 

n 

For this case, T (x(n)) can be shown to be [14], 

T(x(n))= x(n)^H(H^H)r'H^x(n) (15) 
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2 

for the case when o- is assumed to be constant. 

CO 

The optimal estimate 9 for the detected event at n is thus expressed 

as 

e = arg_p(x(n)^(n)) (16) 

e(n) 

5 The IVILE of 9(n)(n) is found by maximizing p(x(n);0(n), or equivalently 

minimizing the IVISE. 

(x(n)-He(n))^(x(n)-He(n))= x(n)'^x(^^ (17) 
Derivation with respect to x(n) yields the optimum 9 , 

9(ii) = (h^h)' H'^x(n) = (h'^h)"* H^x(n) (1 8) 

10 By using the above formulation, it is also possible to derive a 

generalized likelihood ratiodetector based on (15). 
Rate 

The central parameter for classification of cardiac arrhythmias is the 
heart rate. Most arrhythmias are defined in terms of heart rate, although 
15 sometimes with rather fuzzy limits. Consequently, rate should be considered 
in order to improve performance. The RR interval, At, is defined as the 
duration between two consecutive events, 

M,={h,^h,_,)Ts (19) 
where and n,^., denote the occurrence times of the events, and Ts denotes 

20 the sampling period. 

Leader-follower clustering 

The choice of leader-follower clustering is based on a number of 
features which makes it suitable for the present invention, viz. 
o on-line processing (non-iterative), and 
25 o self-learning, i.e., no a priori knowledge of the of clusters is 

needed. 

The starting point is the assumption that an event is present for which it 
should be decided whether it belongs to an already existing cluster or if a new 
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cluster should be initiated. Since, no knowledge is a priori available on which 
rhythms or morphologies to be expected, the chosen algorithm must be self- 
learning. The leader-follower clustering algorithm is constituted by four 
quantities: 

5 1. The event parameter vector, <l>(nk), containing the features of 

the k:th event that occurs at time nk. 

2. The cluster center [}\, and the covariance matrix Zj; that together 
define the i:th cluster. Since both \}\ and li are unknown, they are replaced by 

their estimates ^- and ; respectivey, 
10 3. The metric, dj^ , determine the distance between <t>(nk) and |Ji; 

according to some suitable function. 

4. A rule for adaptation of the cluster parameters for the winning 
cluster. 

Initialization of new clusters 

15 During run-time, a finite number of clusters exists which represent the 

rhythms having appeared until present time. Thus, it is occasionally 
necessary to initialize a new cluster when the existing ones do not sufficiently 
well fit the present event. When the distance function dj^ exceeds a certain 
threshold, ri, it is more likely that the event belongs to a new cluster than to 

20 any of the existing clusters. The selection of n is a tradeoff between cluster 
size and cluster resolution in the sense that choosing a small q will result in 
many clusters with few clustering errors. On the other hand, a large n results 
in few clusters but in more errors. 

The minimum distance between <t>{n) and with respect to both i and 

25 n is 



d^2^ = minc^;(l) i = I = --.fi. +| (20) 

over the search interval K. The corresponding minimum distance 
indices are found as 
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[imin.n^iJ = argmjnd2(|) (21) 
U 

The decision aile based on the comparison of d^j^ and q is 

expressed as 

>r|: Initilialize new cluster 

<r| : Assign to winning cluster (22) 
<pr| : Update winning cluster 



d 2 < 

min 



where 0 < p <1 . When the upper relation in (22) holds, a new cluster is 
initialized by first increasing the number of clusters by one, 1 = 1 + 1, given that 
I < Imax where Uax is the maximum number of clusters, and then initializing the 
new cluster as, 

10 Pi=<I>(nJ (23) 

On the other hand, if I = Imax the algorithm needs to discard one of the 
existing clusters. This can be done by elimination of, e.g., the oldest or the 
"smallest" cluster. By the term "smallest" cluster is meant that cluster which 
most rarely is fitted with a detected cardiac event. 

15 For the middle and lower conditions in (22), the existing minimum 

distance cluster imin is selected as the winning cluster. However, only the 
lower condition results in a cluster parameter update. The reason for 
including such a distinction is that only closely similar events should be used 
for cluster updates in order to reduce contamination. 

20 Mahalanobis distance function 

The distance between the feature vector 0(n) and each cluster , is 
defined as the Mahalanobis distance, which is a normalized Euclidean 
distance in the sense that it projects the parameter vector elements onto 
univariate dimensions by including the inverse covariance matrix Z;'' . Thus, 

25 a feature with a larger variance in 0(n) will be assigned a larger share of the 
hyperspace before normalization compared to that with a lower variance. A 
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consequence of normalization is that the Mahalanobis distance works well on 
correlated data since Ei"* then acts as a decorrelator. 

When searching for the minimum distance, a grid search over n is 
performed. This grid search is necessary since it not only minimizes the 
5 distance but also results in a more accurate fiducial point estimate than what 
would be the case when only considering T(x(nk)). The minimum distance is 

thus found by a grid search with respect to all existing clusters, i = 1 J, and 

all feature vectors within the duration of an event I as defined in (20), 

df(l) = (0(l)-Mi)'t7^(0(l)-p|) (24) 

1 0 Reference feature adaptation 

In order to track changes in the features of the different rhythms, 

adaptation of both p..^ (k) and Z. "} (k) are desirable, here the event index 

imin 

k has been included for clarity. For il^^ (k), an exponentially updated 
average is used: 

15 Mu„(k) = Pu(k-1)+Y£(k) (25) 

where 

£(k) = o(n,J-p;^Jk-l) (26) 

20 The exponential update factor y is confined to the interval 0 < y <^. 

The inverse covariance matrix, ^i^jj(k). is estimated by exponential 

averaging of the new cluster difference matrix £(k)8^(k) using the update 
factor (l-a), 

t : ^ (k) = ((a)i/^in (/c - 1) + (1 - aXky {k)Y (27) 

/ min 

25 

Using the matrix inversion lemma 
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A = B-' +CD-' (28) 
A"' = B - BC(o + C'^BC)"'C^B (29) 

and pairing the terms in (27) with the ones in (28), 

5 A = i,„,„(/f) (30) 

B-'a±.^,„{k-li) (31) 

C = £(/f) (32) 

D-'=(l-a) (33) 

The inverse ^ii^|ii(k-l)in (27) may be computed without any matrix 

1 0 inversions as, 

a-'t rHk- iy{ky{k)a-'t rHk- 1) 

^ A r ^ 'i^ A y V /min^ / \ / \ / /min^ ^ 

(l-ar+eM/tyz^-IC/c-lK*) 

By utilizing the matrix inversion lemma, the computational complexity of 
the operation is reduced from O(P^) to O(P^). 
Algorithm initialization 

15 In leader-follower clustering, clusters are initialized as they become 

needed. This feature is convenient since it does not automatically introduce 
any unused clusters. However, it also puts demands on the algorithm to be 
able to create new clusters when necessary, and also to terminate clusters 
either not used for long or with only a few events. Initially, the total number of 

20 clusters, I, is equal to one. The algorithm is initialized by assigning the 
parameter vector <D(n,), which maximizes the test statistic in (14) of the first 
event, to the first cluster, cf. (23), 



P,=<d(/70 (35) 



25 
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For the general case, ^{n^) is composed of a subset of the 
representative functions together with the preceding RR-interval, 



o(nJ = 



(36) 



5 Whoro where §5(71^) is a subset of the most discriminating elements in 

Not e It should be noted that an event is defined by its depolarization 
wave. Consequently, At is not included in the arrival time estimation of n^. , 
but instead computed afterwards. The time continuous notation of At^ is 
10 preferred since it results in a suitable magnitude similar to the normalized 
morphological information in . 

The inverse correlation matrix estimate is initialized in the same way for 
all clusters; a simple solution is to set it equal to a scaling of the identity matrix 
I. 

15 l;-.l = 5I (37) 

where a is a design parameter. The complete clustering algorithm for 
organized events is presented in the flow chart in Fig. 2. 
Reduced-complexity clustering 

In order to develop a more efficient algorithm in terms of performance 
20 versus power consumption and to evaluate the power consumption itself, a 
simple measure of the computational complexity can be used, namely, the 
total number of multiplications. This quantity represents a much more 
complex operation than do additions. In this algorithm, where most 
operations are of the nature "multiply-accumulate", the number of additions is 
25 of the same order as the number of multiplications and may thus be neglected 
without significant loss of accuracy. 
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In order to reduce the computational complexity, focus is put on 
reducing the number of multiplications. The dominant contributions of 
multiplication operations are found in (24) and (34) which require P(P + 1)and 
P/2 (3P + 5) multiplications, respectively, considering certain symmetry 

5 properties of i^l. Furthermore, one division is required in (34). However, 

according to (20), (24) is performed IK times per event while (34) is performed 
only once per event. 

Based on the above performance figures, a few approximations can be 
identified: 

10 o to use only the peak(s) in T(x(n)) instead of a complete grid 

search, 

o to use a simplified , and 

o to use a likelihood based search sequence over the clusters, 
i.e., to start with the most likely cluster and to stop the search if 
15 a sufficiently small distance is found. 

Simplifying the grid search from spanning both samples and clusters in 
(20) to span over only clusters, 

di^in=«V"^5W (38) 

the number of multiplications may be reduced by a factor K to IP(P +1). 
20 Due to sensitivity in T(x(n)), the feature vectors resulting in the peaks for two 
different events may differ significantly, this simplification is likely to result in 
more clusters. A useful compromise may instead be to use the coefficients 
from, e.g., the 3 largest peaks in T(x(n)) resulting in 3IP(P + 1) multiplications 
per event. 

25 Since a cardiac event lasts for longer time than one sample those 

samples which give the filter coefficients which are most similar to the cluster 
reference are determined. A grid search over 40 msec is preferably made. In 
this way coefficients of greatest importance locally are determined. For this 
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decision T(x(n)) is considered and samples generating a peak are chosen, 
i.e. T(x(n-l)) < T(x(n)) > T(x(n + l)), since they indicate the probability for the 
presence of an event. 

Another simplification, based on the approximation of orthogonality 
5 between the different wavelet scales as well as the RR interval, is to simplify 



where diag(A) returns the diagonal elements of the square matrix A. By using 
this approximation, the number of multiplications used for the estimation of 

10 Z".! is reduced to 3P. Additionally, the distance computation in (20) is 

simplified and may be reduced to 2IKP multiplications per event. 

A reasonable assumption is often that successive events originate from 
the same rhythms. Considering this knowledge, it would be sufficient to 
compute the distances for the previously selected cluster. In doing so, the 
1 5 number of multiplications in (38) are reduced even further. 

Table 1 presents the different detector versions as defined by their 
distinguishing features and shows computational complexity for the different 
versions of the clustering algorithm. 



s -.1 



by only including its diagonal elements in the adaptation, 




(39) 



Features 1 



Version 



i min 



Search alg. 



Complexity C 



Full 



Interval 



IKP(P+ 1)+ -(3P + 5) 



Diagonal 



Interval 



2IKP+ 3P 



Full 



3 peak 




Diagonal 



3 peak 



6IP+3P 
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The total computational complexity, C, as reflected by the number of 
multiplications for the different algorithm versions, is presented in Fig. 3 as a 
function of 1. In figur e Fig. 3 the solid line shows an algorithm version A, 
5 dashed line a version B, dotted line a version C, and a version D by dash- 
dotted line), for a feature vector with (a), P = 4, and (b), P = 7 elements. 
Results 

The results are obtained by studying the performance of the following 
quantities: 

10 o algorithm versions, as presented in Table 1 , and 

o noise tolerance, for noise-free signals and for signals with 

background noise of 50|jV RMS. 
The following parameter settings have been used (unless otherwise 
stated), 

15 

a = 1 .05 y = 0.025 S=50 K=40 (49) 

It is noted that a'^is chosen to offer faster adaption than y • The reason 
for that is that the initial estimate t^,^ is likely to be less accurate than the 

A 

20 initial estimate |Jj . Also, S is chosen to have the same order of magnitude 

as the steady state eigenvalues . 

It should be noted that the different algorithm versions are not fully 
comparable in terms of performance for a specific n due to the differences in 

distance computation. In versions B and D, where a diagonal Z is used, 

25 the lack of non-diagonal information results in a nonorthogonal distance which 
is larger than the orthogonal one. Since versions C and D make use of a 
limited search, the minimum distance found may differ from the global 
minimum distance for the event. Both these algorithmic differences imply an 
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minimum distance for the event. Both these algorithmic differences imply an 
increase in clustering quality for a certain value of q. however, at the expense 
of more introduced clusters. 
Evaluation measures 
5 The two main quantities evaluated are clustering performance and 

computational complexity, see Fig. 3; these two quantities are in general 
conflicting. In the evaluation, a cluster is assigned to that cardiac rhythm 
which contains the most events in the cluster. This rhythm is denoted as the 
dominant rhythm within the cluster. With respect to the dominant rhythm, the 

10 cluster is defined as a correct cluster, whereas it is erroneous to all other 
rhythms. If a rhythm is found to be dominant in more than one cluster, such 
clusters are first merged in the performance evaluation. The number of 
events in the correct cluster which belongs to the i:th dominant rhythm is 
equal to No{\), while the number of events belonging to any other false rhythm 

15 in the cluster is equal to NF(i). The number of events of the dominant rhythm 
which are not classified in a correct cluster, i.e., missing, is equal to NM(i). 

The performance of the algorithm is evaluated in terms of probability of 
a correct clustering of an event, Pcc(i), and probability of a dominant event in 
a cluster, PDE(i). The first parameter is expressed as the share of correctly 

20 clustered events within a rhythm and is, using the above parameters, defined 
by 

Pcc(i) = —J^ (41) 

N^(i)+NJi) 

The second parameter may be expressed as the share of dominant 
rhythm events within a cluster, and is defined as 

25 p^^(/) = _NDfi) (42) 

Np(i)-hNJi) 

For the case when a rhythm completely lacks a correct cluster, PDE(i) is 
undefined; the rhythm is then excluded from subsequent statistical 
computations. Averaging Pcc(i) andPoeCi) over all clusters results in the 
global performance measures Pec and Pde, respectively. 
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It should be pointed out that Pcc(i) and Poe{\) reach their maximal value 
of 1 when r\ is sufficiently small such that the number of clusters equals the 
total number of beats evaluated. This is, of course, a highly undesirable 
solution although performance, as expressed in (41) and (42), will be 
5 excellent. For this reason, the total number of clusters, I, is a crucial 
parameter to be considered. Here, the average I is used together with the 
minimum and maximum number of clusters for a case, Imin. and Uax. 
respectively. 

The power consumption of the algorithm is an important parameter 
10 which determines the pacemaker life span. In this study, power consumption 
is approximated by the computational complexity defined above as the total 
number of multiplications. As shown, the computational complexity depends 
mainly on three parameters: the feature vector length, P, the search interval, 
K and the cluster size, I. 
15 Noise-free signal clustering performance 

F i gur e A pr e s e nts Figs. 4a and 4q illustrate clustering performance in 

terms of Pde and Pec, see Fig. A (a) 4a, and Imin. Tand Uax, see Fig. 4 (b). 
Thus in figur e 4 (a) Fig. 4a clustering performance is shown as depending on 

T for Pec (dark bars) and Pde (bright bars) for noise-free signals. In figur e 
20 4(b) Fig. 4b . the spread between the different cases is shown as, from left to 
right, Imin. I and Lax- The presented algorithm versions are found in Table 1 
and three values of I have been chosen for comparison; 3, 4 and 5. Versions 
A and B perform similarly for all three cases, and achieve Pde = 1 and Pec = 1 

for T = 4 and 1=5, respectively, by creating an acceptable number of 
25 clusters. However, it can be seen from Fig. 4 (b) 4b that a large difference in 
the number of initialized clusters between different cases is present. For 

version B, a slight Increase in both T and Imax is observed for T = 4. However, 
this increase is more due to an unfortunate step-like behaviour in the results 

for the given T than for any significant decrease in performance. The values 
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of n used in Rg-. — 4 Figs. 4a and 4b for the different versions are shown in 
Table 2. 

Table 2: Values of ri used In Fig. 4. 

Algorithm version 



T 


A 


B 


C 


D 


5 


3.6 


4.2 


5.4 


5.8 


4 


4.6 


4.8 


7.0 


7.2 


3 


9.6 


9.6 


12.0 


10.8 



5 

Versions C and D perfomri slightly worse than do versions A and B. 
Contrary to the latter ones, neither version achieves Pde = 1 or Pec = 1 for the 

presented T , see figur e 4 (a) Fig. 4a . 

Figur e 5 (a) pr o s o nts Fig. 5a shows the clustering perfornnance for 
10 version A in detail and its dependence on r|- The most noteworthy result is 
that both Pde = 1 and Pec = 1 for n < 7. It is also clear that clustering 
performance deteriorates rapidly for n >10. The increase in Pec for very high 
values of n is due to that not all rhythms are allocated to a dominant cluster 
and are therefore disregarded in the performance computation. 

15 In Fig. 5 (b) 5b . Imin, \ and Uax are presented as a function of the 

clustering threshold. Thus f i gure 5(a) Fig. 5a shows clustering performance in 
terms of Pde (solid line) and Pec (dashed line) as a function of q for noise-free 

signals, and figure 5(b) the corresponding Imin. T and Imax- For n < 3> f^e 

number of initiated clusters increases rapidly with decreasing ri. Within the 
20 interval 3 < q <7, all rhythms initiate at least one cluster, while for n >7, this is 

not true for all cases. Removing the "worst case", the above is true for the 

other cases for q <10. 

The clustering performance is exemplified for one set of parameters in 

Fig.6 using r\ = 7, Correct clustering with minimal number of clusters is 
25 achieved for two cases while an extra cluster is initialized for two cases. The 
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reason for the extra cluster is that, for case 2, the temporal search interval is 
chosen too small for the third morphology from the left resulting in an extra 
cluster. For case 5 an actual difference in morphologies both on the up and 
down slopes of the dominant peak is discernible between the first and fourth 
5 clusters from the left. 

Clustering for a concatenated electrogram is presented for case 3 in 
Fig. 7, where three distinct rhythm classes can be discerned, viz. normal 
sinus rhythm followed by supraventricular tachycardia and atrial flutter. The 
EGM is shown together with the clusters, represented by o. x and +, 
10 respectively, assigned for each event. The different rhythms result in three 
clusters. 

Figure 8 is a flow chart illustrating the function in broad outline of an 
embodiment of the apparatus according to the invention. At block 2 cardiac 
event features are extracted in the form of wavelet coefficients, and the event 
15 is detected, at 4, 6. At block 8 is checked whether the detected event is 
member of a labelled cluster. If so, the event is added to a class, at 10, and 
actions associated with that class are performed, at 12. 

If the event is not member of a labelled class it is checked if it is a 
member of an existing cluster, at 14. If so, the event is added to a class, at 
20 16, and it is checked if it is possible to label the cluster, 18, and if so the 
cluster is labelled, at 20. 

If the event is not a member of an existing cluster, block 14, a new 
cluster is created as described above fitted to the detected cardiac event, at 
22. 

25 Thus according to the invention clustering events in the EGM is 

performed for use in implantable CRM devices, like heart stimulators. The 
invention is based on feature extraction in the wavelet domain whereupon the 
features are clustered based on the Mahalanobis distance criterion. 
According to advantageous embodiments of the apparatus according to the 

30 invention simplifications of the technique is proposed in order to reduce 
computational complexity to obtain a more implementationally feasible 
solution. 
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If the apparatus according to the invention is to be used for longer 
periods of data analysis, large clusters, although old, may be desirable to be 
kept in some way, while the oldest cluster may be selected to be removed if 
the application is based on shorter time frames. Also, due to the short data 
5 lengths, testing of such algorithms would be of limited value. 

By combining detector/clusterer with labelling rules of a classifier a 
complete detector/classifier is obtained with the possibility to more accurate 
therapy. 

The labelling need not be done in real time and probably more than 
10 one event will be needed to label a cluster. Once the cluster has been 
labelled using clinical terms, the actions associated with the particular class 
will be carried out immediately, i.e. in real time. Thus, the rules needed to 
label the cluster are not used in identifying the event itself. 

The rules used to label the clusters are based on characteristics of the 
15 different possible events. Instead of the exact rules, the characteristics are 
consequently described. 

Figur e Fig. 9 is a block diagram of an embodiment of an implantable 
heart stimulator provided with the apparatus according to the invention. 
Electrodes 30, 32 implanted in the heart 34 of a patient are connected by a 
20 lead 36 to an AID converter 40, via a switch 38 serving as overvoltage 
protection for the A/D converter 40. In the A/D converter 40 the signal is A/D 
converted and the digital signal is supplied to a wavelet detector 42. 

The detector 42 decides whether a cardiac event is present or not as 
described earlier. Wavelet coefficients are calculated as well. Parameters of 
25 the detector 42 are programmable from the stimulator microprocessor 44. At 
the detection the coefficients and the RR information are forwarded to the 
clusterer 46 in which it is determined to which cluster the detected cardiac 
event belongs, as described previously. The clusterer 46 is preferably of a 
leader-follower type and also the cluster parameters are programmable from 
30 the microprocessor 44. 

By the microprocessor 44 suitable therapy is decided depending on 
assigned cluster for the detected event and possible a priori knowledge about 
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arrhythmia associated with the cluster in question. Thus it is possible to 
distinguish e.g. ventricular tachycardia from a sinus tachycardia by comparing 
the parameters with a known normal sinus rhythm. Parameters of the sinus 
tachycardia are then supposed to be similar to those of a normal sinus 
5 rhythm, whereas parameters of ventricular tachycardia differ significantly. 

As an alternative the decision rules can be trained from a number of 
rhythms and the resulting rules are then used on test data, see Weichao Xu et 
al., "New Bayesian Discriminator for Detection of Atrial Tachyarrhythmias", 
DOI:10.1 161/01 .CIR.0000012349.14270.54, pp.1472-1479, January 2002. It 
10 is then possible to decide that a certain position of the cluster indicates e.g. a 
sinus rhythm, etc. Also this technique can be based on analysis of the feature 
vector for a cluster, and it is possible to decide if the beat is broad or narrow, 
large or small, or if the rhythm is regular or irregular. 

The stimulator in figur e Fig. 9 also includes a pulse unit 48 with 
15 associated battery 50 for delivery of stimulation pulses to the patient's heart 
34 depending on the clustering evaluation. 

The implantable stimulator shown in figur e Fig. 9 includes telemetry 
means 52 with antenna 54 for communication with external equipment, like a 
programmer. 

20 Although modifications and changes mav be suggested bv those 

skilled in the art, it is the intention of the inventor to embody within the patent 
warranted hereon all changes and modifications as reasonably and properlv 
come within the scope of his contribution to the art. 
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